function [p,e,t] = new3d_grid(N,limits)
% [p,e,t] = new3d_grid(N,limits)
%
% Build a 3D structured grid using N(1) subdivisions in x, N(2)
% subdivisions in y and N(3) subdivisions in z. The domain is the cube
% [limits(1),limits(2)]x[limits(3),limits(4)]x[limits(5),limits(6)].
% If all the N(:) are even, it is guaranteed that each element has at
% least one internal node.

 nv = (N(1)+1)*(N(2)+1)*(N(3)+1);
 p = zeros(3,nv);
 for i=1:N(1)+1
   for j=1:N(2)+1
     for k=1:N(3)+1
       ijk = (i-1)*(N(2)+1)*(N(3)+1) + (j-1)*(N(3)+1) + k;
       p(:,ijk) = [i,j,k]' - 1;
     end
   end
 end
 for i=1:3
   j = (i-1)*2 + 1;
   k = j+1;
   p(i,:) = ((limits(k)-limits(j))/N(i))*p(i,:) + limits(j);
 end

 ne = 6*N(1)*N(2)*N(3);
 t = zeros(4,ne);
 e = [];
 for i=1:N(1)
   if(2*floor(i/2)==i)
     oj = 'south';
   else
     oj = 'north';
   end
   for j=1:N(2)
     if(2*floor(i/2)==i)
       ok = 'down';
     else
       ok = 'up';
     end
     for k=1:N(3)
       ijk = (i-1)*N(2)*N(3)*6 + (j-1)*N(3)*6 + (k-1)*6;
       iv(1) = (i-1)*(N(2)+1)*(N(3)+1) + (j-1)*(N(3)+1) + k;
       iv(2) = (i-1)*(N(2)+1)*(N(3)+1) + (j-1)*(N(3)+1) + k+1;
       iv(3) = (i-1)*(N(2)+1)*(N(3)+1) + (j-1+1)*(N(3)+1) + k;
       iv(4) = (i-1)*(N(2)+1)*(N(3)+1) + (j-1+1)*(N(3)+1) + k+1;
       iv(5) = (i-1+1)*(N(2)+1)*(N(3)+1) + (j-1)*(N(3)+1) + k;
       iv(6) = (i-1+1)*(N(2)+1)*(N(3)+1) + (j-1)*(N(3)+1) + k+1;
       iv(7) = (i-1+1)*(N(2)+1)*(N(3)+1) + (j-1+1)*(N(3)+1) + k;
       iv(8) = (i-1+1)*(N(2)+1)*(N(3)+1) + (j-1+1)*(N(3)+1) + k+1;
       [tt,ex,ey,ez] = new3d_grid_singlecube(iv,[ok,'_',oj]);
       t(:,ijk+1:ijk+6) = tt;
       if(i==1)
         e = [e,ex(:,:,1)];
       endif
       if(i==N(1))
         e = [e,ex(:,:,2)];
       endif
       if(j==1)
         e = [e,ey(:,:,1)];
       endif
       if(j==N(2))
         e = [e,ey(:,:,2)];
       endif
       if(k==1)
         e = [e,ez(:,:,1)];
       endif
       if(k==N(3))
         e = [e,ez(:,:,2)];
       endif
       if(strcmp(ok,'up')==1)
         ok = 'down';
       else
         ok = 'up';
       endif
     end
     if(strcmp(oj,'north')==1)
       oj = 'south';
     else
       oj = 'north';
     endif
   end
 end

return

